Spectral Graphs & Poly2Graph — A Hands‑On Tutorial

Author

Yiyang Liu

Published

November 10, 2025

1 Introduction

This tutorial is for the ones don’t have physics background, introducing the concepts of spectral graphs and demonstrates how to use the Poly2Graph library to create and manipulate spectral graphs in Python.

To fully benefit from this tutorial, you should have a basic understanding of Python programming and mathematical background on linear algebra and complex functions.

2 Learning objectives

By the end of this tutorial, you will be able to:

  • Explain what a spectral graph is (a graph traced by the energy spectrum of a 1‑D crystal on the complex plane under open boundaries).
  • Understand the key mathematical objects that appear in modern treatments: the characteristic polynomial \(P(z,E)\), the spectral potential \(\Phi(E)\), and the density of states \(\rho(E)\).
  • Describe prior pain points (manual inspection, numerical bottlenecks, fragile extraction) and the need for an automated tool.
  • Use Poly2Graph to go from a symbolic Hamiltonian/characteristic polynomial to its spectral graph, and analyze/visualize the result.

image grid of ChP
Figure 1. Concept map

The Figure 1 shows the concept map of this tutorial, (a) starting from 1-D crystal Hamiltonian \(H(z)\) in momentum space, or its characteristic polynomial \(P(z,E) = z^2 + 1/z^2 + 10Ez - E^4\), (b) calculating the spectral potential \(\Phi(E)\) (Rokin Function) through the roots of \(P(z,E) = 0\), (c) the density of states \(\rho(E)\), obtained as the Laplacian of \(\Phi(E)\) and (d) visualizing the spectral graph through sg.spectral_graph() based on NetWorkX.

3 What is a spectral graph?

Consider a simple 1‑D crystal (a line of repeating “cells” in Figure 1(a)), electrons in a cell are possible hopping to neighbor cells (the “arrow” terms in Figure 1(a)). It is a classic scenario in non-Hermitian system. A good visualized way to uncover the system’s energy structure is to plot a spectral graph. What actually is a spectral graph? Here is a specific explanation.

To seize the concept accurately, we start from a linear algebra viewpoint: studying a family of matrices built on local hoppings (couplings) between neighboring cells.

  1. Non-Hermitian simply means the matrix \(H\) is not equal to its conjugate transpose \(H^\dagger\). This can happen if forward and backward hoppings differ or if we have gain/loss factors (e.g. non-reciprocal hoppings or some sort of on-site potentials). Consequently its eigenvalues need not be real, they live in \(\mathbb{C}\).

  2. Two boundary conditions cause two global matrix structures:

    • Periodic Boundary Conditions (PBC): allowing hoppings between boundaries, the matrix “wraps around” (i.e., the first and last cells are coupled) — it is (block) circulant.
    • Open Boundary Conditions (OBC): turning off boundary hoppings, the wrap is cut — the matrix is (block) banded Toeplitz.

  3. Under PBC (circulant case) eigenvalues are easy to understand: a circulant (or block circulant) matrix is diagonalized by Discrete Fourier Transform (DFT) from real space to momentum space. Writing a small “Bloch” matrix \(h(z)\) with \(z=e^{ik}\), sampling \(k\in[0,2\pi)\) makes the eigenvalues trace smooth closed curves or loops in the complex plane.

  4. Under OBC (Toeplitz case) the same local rules yield a different global matrix. In the thermodynamic limit (number of cells \(N\to\infty\)), the OBC spectrum generally does not fill those PBC loops. Instead, it is supported on thin curves (arcs and sometimes closed loops) in the complex \(E\)‑plane. Connecting these curves yields a planar network: the spectral graph.

In conclusion: a spectral graph is the OBC spectrum in the thermodynamic limit (\(N\to\infty\)) viewed as a network of curves in the complex \(E\)‑plane.

For more mathematical details, see Section 7.

4 Why interesting?

Now we have a basic understanding of what a spectral graph is. However, which domains can benefit from spectral graphs? Why are spectral graphs interesting? Here are some views:

Within Non‑Hermitian systems:

Non-Hermitian spectra can realize a plethora of graph connectivities than standard topological labels (like \(\mathbb{Z},\mathbb{Z}_2\)): stars, flowers, braids, etc., with novel transitions that change the number of endpoints/junctions/loops (no Hermitian analog). This is naturally expressed by writing dispersions as bivariate Laurent polynomials \(P(E,z)\) and studying how their imaginary‑momentum (“inverse skin depth” \(\log|z|\)) surfaces intersect across the \(E\)‑plane.

Beyond Non‑Hermitian (algebra‑to‑graph):

Spectral graph is a general algebra-to-graph bridge. Because the equation \(\Phi(E) = - \log|a_q(E)| - \sum_{i=p+1}^{p+q}\log|\beta_i(E)|\) uses only \(P(E,z)\), thus

  1. Any bivariate Laurent polynomial \(P(z,E)\) has a spectralgraph.
  2. Any univariate polynomial \(P(z)\) can be viewed as one-band Bloch Hamiltonian; and any vector can be treated as a symmetrised coefficient list of a univariate polynomial.
  3. Any matrix can be decomposed into a product of one-band Hamiltonian matrices, which can form a multi-band Bloch Hamiltonian, so in general it has a multiset of spectral graphs.

You can get a universal topological fingerprint for many algebraic objects—making graph learning a front-door to problems in linear algebra, signal processing, photonics, circuits, and more.

5 A quick start: One‑band example

We’ll start from a simple one‑band model and obtain its spectral graph.

Characteristic polynomial

The model’s characteristic polynomial can be written as \(P(z,E) = z^4 - z^2 - z^{-1} - E\) (This corresponds to a scalar Bloch Hamiltonian \(h(z)=z^4-z^2-z^{-1}\), with \(z=e^{ik}\)).

# If running on devices without essential libraries, uncomment the next line:
# !pip install -U poly2graph sympy networkx matplotlib

import sympy as sp
import numpy as np
np.set_printoptions(threshold=20)
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import networkx as nx

# Poly2Graph (must be installed in your environment)
# If this import fails, run the pip command above.
import poly2graph as p2g

Please ensure your version of Python >= 3.11 , you can check the installation:

print("Poly2Graph version:", p2g.__version__)
Poly2Graph version: 0.2.0
# Always remember to define common symbols first
k = sp.symbols('k', real=True)
z, E = sp.symbols('z E', complex=True)

The valid input formats to initialize a p2g.SpectralGraph object are:

  1. Characteristic polynomial in terms of z and E:

    • as a string of the Poly in terms of z and E
    • as a sympy’s Poly (sympy.polys.polytools.Poly) with {z, 1/z, E} as generators

  2. Bloch Hamiltonian in terms of k or z

    • as a sympy Matrix in terms of k
    • as a sympy Matrix in terms of z

All the following characteristics are valid and will initialize to the same characteristic polynomial and therefore produce the same spectral graph

char_poly_str = '-z**-1 - E - z**2 + z**4'

char_poly_Poly = sp.Poly(
    -z**-1 - E - z**2 + z**4,
    z, 1/z, E # generators are z, 1/z, E
)

phase_k = sp.exp(sp.I*k)
char_hamil_k = sp.Matrix([-phase_k**-1 - phase_k**2 + phase_k**4])

char_hamil_z = sp.Matrix([-z**-1 - E - z**2 + z**4])

Here we use the string format as the initialization example.

sg = p2g.SpectralGraph(char_poly_str, k=k, z=z, E=E)
# print the number of energy bands
print('Bands:', sg.num_bands)

# print the max hopping range (p,q)
print('Max hop (right):', sg.poly_p, 'Max hop (left):', sg.poly_q)

# print the [Ereal_min,Ereal_max,Eimag_min,Eimag_max] window
print('Spectral window:', 
      dict(zip(['Re(E)_min', 'Re(E)_max', 'Im(E)_min', 'Im(E)_max'], 
               map(float, sg.spectral_square))))
Bands: 1
Max hop (right): 1 Max hop (left): 4
Spectral window: {'Re(E)_min': -2.3113904784272243, 'Re(E)_max': 1.6398687049402647, 'Im(E)_min': -1.9756295916837434, 'Im(E)_max': 1.9756295916837456}

We will see Poly2Graph automatically compute set of properties next:

sg.ChP

\(\displaystyle \operatorname{Poly}{\left( z^{4} - z^{2} - \frac{1}{z} - E, z, \frac{1}{z}, E, domain=\mathbb{Z} \right)}\)

Bloch Hamiltonian:

  • For one-band model(The exponent of \(E\) is 1), it is a unique, \(1 \times 1\) matrix (scalar).
sg.h_k

\(\displaystyle \left[\begin{matrix}e^{4 i k} - e^{2 i k} - e^{- i k}\end{matrix}\right]\)

  • or in \(z\) format(recall that \(z = e^{ik}\)).
sg.h_z

\(\displaystyle \left[\begin{matrix}- \frac{- z^{5} + z^{3} + 1}{z}\end{matrix}\right]\)

Frobenius Companion Matrix of \(P(z, E)\):

  • The Frobenius companion matrix is a special matrix whose eigenvalues are exactly the roots of the characteristic polynomial \(P(z, E)\) for fixed \(E\).
  • In Poly2Graph, sg.companion_E gives this matrix (with \(E\) as a parameter, \(z\) as variable).
  • Purpose:
    • Efficiently computes all \(z\)-roots for any \(E\) using linear algebra (eigenvalues), instead of slow generic root-finding algorithms.
    • These \(z\)-roots are essential for calculating the spectral graph, GBZ, spectral potential, and density of states.
sg.companion_E

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 1\\1 & 0 & 0 & 0 & E\\0 & 1 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 1\\0 & 0 & 0 & 1 & 0\end{matrix}\right]\)

It is a key tool for automating and accelerating spectral graph extraction in Poly2Graph.

5.1 Spectral potential \(\Phi(E)\), density of states \(\rho(E)\), and skeleton

We compute \(\Phi(E)\) and \(\rho(E)\) on a grid of complex energies, then binarize \(\rho(E)\) to get the skeleton of the spectral graph. The phi and dos attributes store \(\Phi(E)\) and \(\rho(E)\) values on the grid, while skel is a boolean array indicating where \(\rho(E)\) has supported.

phi, dos, skel = sg.spectral_images()

fig, axes = plt.subplots(1, 3, figsize=(8, 3), sharex=True, sharey=True)

# 1) Spectral potential
axes[0].imshow(phi, extent=sg.spectral_square, origin='lower', cmap='terrain')
axes[0].set(xlabel='Re(E)', ylabel='Im(E)', title='Spectral Potential')

# 2) Density of States
pmin, pmax = np.percentile(dos, (0.1, 99.9))
# Clip extreme DOS to increase visibility
norm = colors.Normalize(vmin=pmin, vmax=pmax)
axes[1].imshow(dos, extent=sg.spectral_square, cmap='viridis', norm=norm)
axes[1].set(xlabel='Re(E)', title='Density of States')

# 3) Binarized skeleton
axes[2].imshow(skel, extent=sg.spectral_square, cmap='gray')
axes[2].set(xlabel='Re(E)', title='Graph Skeleton')
plt.tight_layout()
plt.savefig("assets/1-band-example_si.png")
plt.show()

6 What was hard before (and why we need Poly2Graph)?

Prior to Poly2Graph, extracting spectral graphs was manual and fragile:

  • You had to manually plot long chain spectra, and visually infer the continum graph; for some scenarios of extremely complex spectra, they are plotted slolwly and fragile.
  • Computing the GBZ roots \(\beta_i(E)\) for every sampled energy \(E\) is a bottleneck; naive solvers are slow and memory‑hungry.
  • Even with a DOS image, turning that into a clean graph with correct endpoints, junctions, and loops requires careful morphology (thresholding, skeletonization, node/edge parsing).

Poly2Graph lifts these barriers with an automated end‑to‑end, high‑performance pipeline:

  1. Accepts either \(H(z)\) (as a sympy.Matrix) or \(P(z,E)\) (string or sympy.Poly).
  2. Computes the spectral potential \(\Phi(E)\) from \(z\)‑roots (fast eigen‑solvers on Frobenius companion matrices; optional GPU/TPU acceleration via TensorFlow/PyTorch).
  3. Takes the Laplacian on Rokin function to get DOS, then does adaptive, two‑stage refinement (coarse mask → high‑res on just ~1–5% of the region).
  4. In rare low‑DOS or long‑range multi‑band settings, local numerical instabilities can arise, however, Poly2Graph mitigates them with node‑merging and short‑edge contraction.

This automation has enabled HSG‑12M (11.6M static + 5.1M temporal Hamiltonian spectral multigraphs), demonstrating scale and diversity previously out of reach (Learn more from HSG-12M: A Large-Scale Spatial Multigraph Dataset if interested).

7 Mathematics and Physics

In this section, we summarize the key mathematical objects that appear in modern treatments of spectral graphs: the characteristic polynomial \(P(z,E)\), the spectral potential \(\Phi(E)\), and the density of states \(\rho(E)\).

Algebraic formulation:

Let \(H(z)\) be the Bloch Hamiltonian with \(z=e^{ik}\) and \[ P(z,E)=\det\!\big[H(z)-EI\big]=\sum_{n=-p}^{q} a_n(E)\,z^n . \] For a fixed \(E\), order the \(z\)‑roots of \(P(z,E)=0\) by magnitude \(|\beta_1(E)|\le\cdots\le|\beta_{p+q}(E)|\). The generalized Brillouin zone (GBZ) is characterized by the condition: \[ |\beta_p(E)|=|\beta_{p+1}(E)|, \]

and the corresponding loci of \(E\) form the OBC spectral graph. In particular, the density of states (DOS) in the thermodynamic limit (Number of cells \(N\to\infty\)) is supported only along these loci in the distributional sense.

Potential‑theoretic view:

A convenient “spectral potential” (Ronkin function) can be defined from the \(z\)‑roots; one common choice is: \[ \Phi(E)=-\log|a_q(E)|-\sum_{i=p+1}^{p+q}\log|\beta_i(E)|. \] Where \(a_q(E)\) is the leading coefficient of the characteristic polynomial \(P(z,E)\). Its Laplacian yields the density of states(DOS), \[ \rho(E)=-\frac{1}{2\pi}\,\nabla^2\Phi(E), \] where the “\(-\)” sign depends on the convention used for \(\Phi\). In the distributional sense, \(\rho(E)\) is supported only on the loci of spectral graph. Equivalently, the spectral graph coincides with the ridges of \(\Phi(E)\).

Intuition: treat each eigenvalue \(\epsilon_n\) as a tiny electric charge \(\rho(E) = \lim_{N\to\infty}\sum_n \frac{1}{N} \delta(E-\epsilon_n)\); while \(\Phi(E)\) is an electrostatic‑like potential whose “ridges” trace the spectral graph; \(\rho(E)\) is nonzero only on the graph. Please see the Section 12 below for a concrete illustration.

8 Extract the spectral graph

Continue with the one-band model in the quick start part, we visualize these results by an interactive 3D surface plot of \(\Phi(E)\). Feel free to rotate/zoom the plot to examine the surface and skeleton from different angles.

import plotly.graph_objects as go

E_re = np.linspace(*sg.spectral_square[:2], phi.shape[0])
E_im = np.linspace(*sg.spectral_square[2:], phi.shape[1])
fig = go.Figure(data=[go.Surface(z=phi, x=E_re, y=E_im, 
            opacity=0.6, colorscale='Spectral_r')])
fig.update_layout(
scene=dict(
    aspectratio=dict(x=1.5, y=1, z=.8),
    xaxis_title='Re(E)',yaxis_title='Im(E)',zaxis_title='Phi(E)',
),
title='3D Surface Plot of Spectral Potential',
width=700, height=600
)
fig.show()
Figure 1

sg.spectral_graph() returns a NetworkX MultiGraph with node positions and polyline geometry along each edge.

# Draw the multi-band spectral graph
G = sg.spectral_graph()
pos2 = nx.get_node_attributes(G, 'pos')
plt.figure()
nx.draw_networkx_nodes(G, pos2, node_size=20)
nx.draw_networkx_edges(G, pos2)
plt.axis('equal'); plt.title('Spectral Graph (1-band)')
plt.show()

print(f'#nodes: {G.number_of_nodes()}  #edges: {G.number_of_edges()}')

# Show attributes of a node and an edge
any_node = list(G.nodes(data=True))[0]
print('A node has attributes:', list(any_node[1].keys()))
any_edge = list(G.edges(keys=True, data=True))[0]
print('An edge has attributes:', list(any_edge[3].keys()))

#nodes: 7  #edges: 6
A node has attributes: ['pos', 'dos', 'potential']
An edge has attributes: ['weight', 'pts', 'avg_dos', 'avg_potential']

As you can see above, the spectral graph is a networkx.MultiGraph object.

  • Node Attributes
    1. pos : (2,)-numpy array
      • the position of the node \((\text{Re}(E), \text{Im}(E))\)
    2. dos : float
      • the density of states at the node
    3. potential : float
      • the spectral potential at the node
  • Edge Attributes
    1. weight : float
      • the weight of the edge, which is the length of the edge in the complex energy plane
    2. pts : (w, 2)-numpy array
      • the positions of the points constituting the edge, where w is the number of points along the edge, i.e., the length of the edge, equals weight
    3. avg_dos : float
      • the average density of states along the edge
    4. avg_potential : float
      • the average spectral potential along the edge
node_attr = dict(G.nodes(data=True))
edge_attr = list(G.edges(data=True))
print('The attributes of the first node\n', node_attr[0], '\n')
print('The attributes of the first edge\n', edge_attr[0][-1], '\n')
The attributes of the first node
 {'pos': array([0.88936098, 1.90810319]), 'dos': np.float32(0.25110298), 'potential': np.float32(-0.4699422)} 

The attributes of the first edge
 {'weight': np.float64(2.0950480742544872), 'pts': array([[0.88550233, 1.89652723],
       [0.88550233, 1.89266858],
       [0.88550233, 1.88880993],
       ...,
       [0.81990525, 0.01350528],
       [0.8237639 , 0.00964663],
       [0.8237639 , 0.00192933]]), 'avg_dos': np.float32(0.079228334), 'avg_potential': np.float32(-0.2639912)}

9 A multi‑band example

Now consider a four‑band chain, its polynomial can be written as:

\[ P(z,E) = z^2 + z^{-2} + E z - E^4 \]

A possible Bloch Hamiltonian realization is: \[\textbf{h}(z)=\begin{pmatrix} 0 & 0 & 0 & z^2 + 1/z^2 \\ 1 & 0 & 0 & z \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{pmatrix}\] Poly2Graph can work directly from the string/sympy.Poly as shown below.

# We start with a polynomial string
sg_multi = p2g.SpectralGraph("z**2 + 1/z**2 + E*z - E**4", k, z, E)
# print the number of energy bands
print('Bands:', sg_multi.num_bands)

# print the max hopping range (p,q)
print('Max hop (right):', sg_multi.poly_p, 'Max hop (left):', sg_multi.poly_q)

# print the [Ereal_min,Ereal_max,Eimag_min,Eimag_max] window
print('Spectral window:', 
      dict(zip(['Re(E)_min', 'Re(E)_max', 'Im(E)_min', 'Im(E)_max'], 
               map(float, sg_multi.spectral_square))))
Bands: 4
Max hop (right): 2 Max hop (left): 2
Spectral window: {'Re(E)_min': -1.4055954628228848, 'Re(E)_max': 1.4055954628228862, 'Im(E)_min': -1.4055954628228862, 'Im(E)_max': 1.4055954628228848}
sg_multi.ChP

\(\displaystyle \operatorname{Poly}{\left( z^{2} + zE + \frac{1}{z^{2}} - E^{4}, z, \frac{1}{z}, E, domain=\mathbb{Z} \right)}\)

sg_multi.h_k

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & e^{2 i k} + e^{- 2 i k}\\1 & 0 & 0 & e^{i k}\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

sg_multi.h_z

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & z^{2} + \frac{1}{z^{2}}\\1 & 0 & 0 & z\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

sg_multi.companion_E

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & -1\\1 & 0 & 0 & 0\\0 & 1 & 0 & E^{4}\\0 & 0 & 1 & - E\end{matrix}\right]\)

9.1 Spectral potential \(\Phi(E)\), density of states \(\rho(E)\), and skeleton

Like the one-band example, we also compute \(\Phi(E)\) and \(\rho(E)\) on a grid of complex energies, then binarize \(\rho(E)\) to get a skeleton.

phi2, dos2, skel2 = sg_multi.spectral_images()

fig, axes = plt.subplots(1, 3, figsize=(8, 3), sharex=True, sharey=True)

# 1) Spectral potential
axes[0].imshow(phi2, extent=sg.spectral_square, origin='lower', cmap='terrain')
axes[0].set(xlabel='Re(E)', ylabel='Im(E)', title='Spectral Potential')

# 2) Density of States
pmin, pmax = np.percentile(dos, (0.1, 99.9))
# Clip extreme DOS to increase visibility
norm = colors.Normalize(vmin=pmin, vmax=pmax)
axes[1].imshow(dos2, extent=sg.spectral_square, cmap='viridis', norm=norm)
axes[1].set(xlabel='Re(E)', title='Density of States')

# 3) Binarized skeleton
axes[2].imshow(skel2, extent=sg.spectral_square, cmap='gray')
axes[2].set(xlabel='Re(E)', title='Graph Skeleton')
plt.tight_layout()
plt.savefig("assets/multi-band-example_si.png")
plt.show()

9.2 Extract the spectral graph

# Draw the multi-band spectral graph
G2 = sg_multi.spectral_graph()
pos2 = nx.get_node_attributes(G2, 'pos')
plt.figure()
nx.draw_networkx_nodes(G2, pos2, node_size=20)
nx.draw_networkx_edges(G2, pos2)
plt.axis('equal'); plt.title('Spectral Graph (multi‑band)')
plt.show()

print(f'#nodes: {G2.number_of_nodes()}  #edges: {G2.number_of_edges()}')

# Show attributes of a node and an edge
any_node = list(G2.nodes(data=True))[0]
print('A node has attributes:', list(any_node[1].keys()))
any_edge = list(G2.edges(keys=True, data=True))[0]
print('An edge has attributes:', list(any_edge[3].keys()))

#nodes: 21  #edges: 18
A node has attributes: ['pos', 'dos', 'potential']
An edge has attributes: ['weight', 'pts', 'avg_dos', 'avg_potential']

10 Characteristic Polynomial Class (p2g.CharPolyClass)

For scanning families of polynomials (e.g., varying two complex coefficients across a grid), use p2g.CharPolyClass. This class parallelizes evaluation and is designed to generate many graphs efficiently.

To be specific, p2g.CharPolyClass is:

  • A class of parametrized non-Hermitian lattices
  • Generate spectral properties, including spectral graph in parallel
  • Optimized for computational efficiency and numerical stability

Let us add two parameters {a,b} to the aforementioned multi-band example and construct a p2g.CharPolyClass object.

a, b = sp.symbols('a b', real=True)

cp = p2g.CharPolyClass(
    "z**2 + a/z**2 + b*E*z - E**4", 
    k=k, z=z, E=E,
    params={a, b}, # pass parameters as a set
)
Derived Bloch Hamiltonian `h_z` with 4 bands.

Just like the examples above, we can also check the automatically computed properties:

cp.ChP

\(\displaystyle \operatorname{Poly}{\left( z^{2} + b zE + a \frac{1}{z^{2}} - E^{4}, z, \frac{1}{z}, E, domain=\mathbb{Z}\left[a, b\right] \right)}\)

cp.h_k

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & a e^{- 2 i k} + e^{2 i k}\\1 & 0 & 0 & b e^{i k}\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

cp.h_z

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & \frac{a}{z^{2}} + z^{2}\\1 & 0 & 0 & b z\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

cp.companion_E

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & - a\\1 & 0 & 0 & 0\\0 & 1 & 0 & E^{4}\\0 & 0 & 1 & - E b\end{matrix}\right]\)

To get an array of spectral images or spectral graphs, we first prepare the values of the parameters {a,b}

a_array = np.linspace(-2, 2, 6)
b_array = np.linspace(-1, 1, 6)
a_grid, b_grid = np.meshgrid(a_array, b_array)
param_dict = {a: a_grid, b: b_grid}
print('a_grid shape:', a_grid.shape,
    '\nb_grid shape:', b_grid.shape)
a_grid shape: (6, 6) 
b_grid shape: (6, 6)

You can check your parameter arrays.

a_grid, b_grid
(array([[-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ]]),
 array([[-1. , -1. , -1. , -1. , -1. , -1. ],
        [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6],
        [-0.2, -0.2, -0.2, -0.2, -0.2, -0.2],
        [ 0.2,  0.2,  0.2,  0.2,  0.2,  0.2],
        [ 0.6,  0.6,  0.6,  0.6,  0.6,  0.6],
        [ 1. ,  1. ,  1. ,  1. ,  1. ,  1. ]]))

Note that the value array of the parameters should have the same shape, which is also the shape of the output array of spectral images.

phi_arr, dos_arr, binarized_dos_arr, spectral_square = \
    cp.spectral_images(param_dict=param_dict, device='/cpu:0')
print('phi_arr shape:', phi_arr.shape,
    '\ndos_arr shape:', dos_arr.shape,
    '\nbinarized_dos_arr shape:', binarized_dos_arr.shape)
phi_arr shape: (6, 6, 1024, 1024) 
dos_arr shape: (6, 6, 1024, 1024) 
binarized_dos_arr shape: (6, 6, 1024, 1024)

phi_arr, dos_arr, and binarized_dos_arr are arrays computed by cp.spectral_images(), representing different spectral quantities:

  • phi_arr: The array of spectral potential values, i.e., \(\Phi(E)\), for each parameter point ((a,b), here \(6*6\)) and energy grid (each small sampling grid on the complex \(E\) plain, here \(1024*1024\)). It describes the “potential landscape” over the complex energy plane, with its ridges tracing the spectral graph.
  • dos_arr: The array of density of states values, i.e., \(\rho(E)\), for each parameter point (here \(6*6\)) and energy grid (here \(1024*1024\)). Mathematiccally the laplacian of \(\Phi(E)\), intuitively indicating the positions of the ridges in phi_arr.
  • binarized_dos_arr: The binarized density of states array, obtained by thresholding \(\rho(E)\), which extracts the skeleton (the graph structure) for further analysis.

These arrays typically have shapes like (number of parameter points, ..., H, W), where \(H\) and \(W\) are the height and width of the energy grid. Each parameter point corresponds to one spectral potential/density/skeleton image.

from mpl_toolkits.axes_grid1 import ImageGrid

fig = plt.figure(figsize=(13, 13))
grid = ImageGrid(fig, 111, nrows_ncols=(6, 6), axes_pad=0, 
                 label_mode='L', share_all=True)

for ax, (i, j) in zip(grid, [(i, j) for i in range(6) for j in range(6)]):
    ax.imshow(phi_arr[i, j], extent=spectral_square[i, j], cmap='Spectral_r')
    ax.set(xlabel='Re(E)', ylabel='Im(E)')
    ax.text(
        0.03, 0.97, f'a = {a_array[i]:.2f}, b = {b_array[j]:.2f}',
        ha='left', va='top', transform=ax.transAxes,
        fontsize=10, color='tab:red',
        bbox=dict(alpha=0.8, facecolor='white')
    )

plt.savefig("assets/ChP-example_si.png")
plt.show()

12 SSH example

Consider a SSH (Su–Schrieffer–Heeger) chain, which is an extremely important model for studying the fundamental properties of one-dimensional topological insulators.

SSH connection illustration
Figure 3. SSH chain with alternating hopping strengths

Here solid and dashed lines in the figure mean different hopping strengths. To simplify the problem, the hopping range is limited to one.

The Bloch Hamiltonian reads \(H(z)= \begin{pmatrix} 0 & t_1+t_2 z \\ t_1+t_2 z^{-1} & 0 \end{pmatrix}, \qquad z=e^{ik}.\) In real space this corresponds to a bipartite lattice where intracell and intercell hoppings alternate along the chain.

Using Schrödinger Equation \(H\psi = E\psi\), transposing to \((H - EI)\psi = 0\). Since \(\psi \neq 0\), we take the determinant yields the characteristic polynomial:

\[\begin{aligned} P(z,E) &= \det[H(z)-EI] \\ &=E^{2}-(t_1+t_2 z)(t_1+t_2 z^{-1}) \\ &=E^{2}-(t_1^2+t_2^2+t_1 t_2(z+z^{-1})). \end{aligned}\]

We can obtain the \(z(E)\) by simply solving the roots of \(P(z,E) = 0\): \(z_{\pm}(E) =\frac{\,E^2 - t_1^2 - t_2^2 \,\pm\, \sqrt{(E^2 - t_1^2 - t_2^2)^2 - 4\,t_1^2 t_2^2}\,} {\,2 t_1 t_2\,}.\) here \(z(E)\) and \(\beta (E)\) are the same with different mathematical symbols.

Due to the maximum hopping range to left and right here are both one, we notice that \(|\beta_p(E)|=|\beta_{p+1}(E)|\) present as \(|z_{-}(E)| = |z_{+}(E)|\), thus we can define the GBZ of this SSH model.

Using the equations mentioned above Section 7 , we can also try to obtain the \(\Phi(E)\) and \(\rho(E)\) ourselves.

You can check your answer with mine:

\(\Phi(E) = - \log|t_{1}t_{2}| - \log|z_{+}(E)|\), \(\rho(E) = \frac{1}{2\pi}\nabla^2 \log|z_{+}(E)|.\)

Next, we can use Poly2Graph to check the properties of SSH model and visualize the spectral graph of it. If you are confused about the code, please refer to the previous examples Section 5.

# Always remember to define common symbols first
k = sp.symbols('k', real=True)
z, E = sp.symbols('z E', complex=True)

Here we use the string format as the initialization example.

# here t1, t2 = 1, 1, using P(z,E) = E^2 - (1 + z)(-1 + z^-1)
char_poly_str_ssh = 'E**2 + z - z**-1'
sg_ssh = p2g.SpectralGraph(char_poly_str_ssh, k=k, z=z, E=E)
# print the number of energy bands
print('Bands:', sg_ssh.num_bands)

# print the max hopping range (p,q)
print('Max hop (right):', sg_ssh.poly_p, 'Max hop (left):', sg_ssh.poly_q)

# print the [Ereal_min,Ereal_max,Eimag_min,Eimag_max] window
print('Spectral window:', 
      dict(zip(['Re(E)_min', 'Re(E)_max', 'Im(E)_min', 'Im(E)_max'], 
               map(float, sg_ssh.spectral_square))))
Bands: 2
Max hop (right): 1 Max hop (left): 1
Spectral window: {'Re(E)_min': -1.0484584139607451, 'Re(E)_max': 1.048458413960742, 'Im(E)_min': -1.0484584139607425, 'Im(E)_max': 1.0484584139607447}
sg_ssh.ChP

\(\displaystyle \operatorname{Poly}{\left( z - \frac{1}{z} + E^{2}, z, \frac{1}{z}, E, domain=\mathbb{Z} \right)}\)

sg_ssh.h_z

\(\displaystyle \left[\begin{matrix}0 & - z + \frac{1}{z}\\1 & 0\end{matrix}\right]\)

sg_ssh.companion_E

\(\displaystyle \left[\begin{matrix}0 & 1\\1 & - E^{2}\end{matrix}\right]\)

phi, dos, skel = sg_ssh.spectral_images()

fig, axes = plt.subplots(1, 3, figsize=(8, 3), sharex=True, sharey=True)

# 1) Spectral potential
axes[0].imshow(phi, extent=sg.spectral_square, origin='lower', cmap='terrain')
axes[0].set(xlabel='Re(E)', ylabel='Im(E)', title='Spectral Potential')

# 2) Density of States
pmin, pmax = np.percentile(dos, (0.1, 99.9))
# Clip extreme DOS to increase visibility
norm = colors.Normalize(vmin=pmin, vmax=pmax)
axes[1].imshow(dos, extent=sg.spectral_square, cmap='viridis', norm=norm)
axes[1].set(xlabel='Re(E)', title='Density of States')

# 3) Binarized skeleton
axes[2].imshow(skel, extent=sg.spectral_square, cmap='gray')
axes[2].set(xlabel='Re(E)', title='Graph Skeleton')
plt.tight_layout()
plt.savefig("assets/ssh-example_si.png")
plt.show()